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Abstract — We consider a broad class of interference coordina- 
tion and resource allocation problems for wireless links where 
the goal is to maximize the sum of functions of individual 
link rates. Such problems arise in the context of, for example, 
fractional frequency reuse (FFR) for macro-cellular networks and 
dynamic interference management in femtocells. The resulting 
optimization problems are typically hard to solve optimally even 
using centralized algorithms but are an essential computational 
step in implementing rate-fair and queue stabilizing scheduling 
policies in wireless networks. We consider a belief propagation 
framework to solve such problems approximately. In particular, 
we construct approximations to the belief propagation iterations 
to obtain computationally simple and distributed algorithms with 
low communication overhead. Notably, our methods are very 
general and apply to, for example, the optimization of transmit 
powers, transmit beamforming vectors, and sub-band allocation 
to maximize the above objective. Numerical results for femtocell 
deployments demonstrate that such algorithms compute a very 
good operating point in typically just a couple of iterations. 

Index Terms — Interference coordination, cellular systems, 
wireless communications, belief propagation, femtocells. 



I. Introduction 

Interference coordination has re-emerged as a fundamen- 
tal challenge for next-generation cellular wireless systems. 
Traditional macrocellular deployments are likely to be sup- 
plemented with smaller femtocells and relays, with mixtures 
of restricted and open access, often deployed in an ad hoc 
manner QJ, (2). Such deployments may create much stronger 
and highly variable (in time and space) interference conditions 
than those experienced in current macrocellular networks, and 
traditional cellular power and rate control may not be adequate 
||3). To address this challenge, a key focus of the current 
3GPP LTE-Advanced standardization efforts is on the design 
of an interference coordination framework for such unplanned 
cellular deployments of base-stations with widely different 
transmission powers [4|. Current release of the LTE specifi- 
cation [5 1 provides simple methods for inter-cell interference 
coordination (ICIC), see e.g., [6|. Along with the design of 
mechanisms for ICIC, algorithms to exploit such mechanisms 
are an active area of research as well, see for example, J7J, (SJ. 

Mathematically, interference coordination is a complex dis- 
tributed optimization problem involving scheduling decisions 
at the transmitters of multiple interfering links. In this work, 
we consider a general linear mixing interference model where 
the scheduling decisions in each link are represented as vector 
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(e.g., transmission powers on different sub-bands in frequency, 
beamforming weights) and the interference on each link is 
a linear combination of the scheduling vectors on the other 
links. Associated with each link at a given time is a utility 
function which describes the benefit to a link as a function 
of the scheduling vector from the serving transmitter and 
interference from the other transmitters. The linear mixing 
model is extremely general and can apply to a large class 
of interference models and objectives. Computing maximum 
weighted matching for queue stability |9| and maximization 
of sum utility of average rates for fairness (jTOj are special 
cases of this formulation. 

In the past few years, algorithms for special cases of 
the above problem have been extensively studied. In many 
cases, algorithms with provable desired properties have been 
obtained - for example, there is a rich literature on distributed 
power control methods to achieve a desired SINR for each 
link pT| and to maximize a certain class of utility functions 
of SINRs [ 12 1, approximation algorithms for maximum weight 
matching for combinatorial interference model were obtained 
in |T3j, fl4) , efficient methods to compute optimal beam- 
forming vector for multiuser downlink fT3) and maximizing 
sum utility of SINRs on uplink |16|, stabilizing policies for 
collision sense multiple access (CSMA) type of models based 
on simulated annealing were derived in p7| . More generally, 
heuristic algorithms have been constructed to solve certain 
specific problems approximately in, for example, ]T8) , [19|. 
While these algorithms perform well in practice in spite of no 
provable guarantees, the insights and approximations used to 
obtain these algorithms are very specific to the problem under 
consideration. 

In this paper, we make the following contributions: 
BP Framework: We consider a belief propagation (BP) frame- 
work for a very general wireless scheduling and interference 
coordination optimization problem. The underlying optimiza- 
tion problem is posed as a problem of estimating marginals of 
a joint probability distribution. BP provides a systematic and 
general approach to obtain distributed algorithms; it can be 
used with arbitrary nonlinear utility functions and scheduling 
vectors sets, which enable the algorithm to be applied a range 
of complex scheduling problems including power control, 
subband scheduling and distributed beamforming. Also, while 
we do not obtain any theoretical guarantees, in practice a 
few iterations of BP generate a good operating point. Thus, it 
typically has faster convergence than gradient based algorithms 
(e.g., 10s of iterations in [16]) or simulated annealing [20]. 
Approximation Algorithms: It is well known that implementing 
BP for distributed optimization problems entails high compu- 
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tational complexity and communication overhead. Exploiting 
the linear mixing interference model and applying Gaussian 
and first-order approximations similar to (2"T)-j2"5), we de- 
velop an approximate BP method that has low complexity, 
distributed implementation and minimal messaging. Along 
many links, messages can be carried in small payloads and 
can be broadcast without separate unicast transmissions, which 
is particularly crucial for wireless systems. Moreover, the 
resulting algorithm has a natural interpretation has a "soft" 
RTS / CTS type handshaking. The approximate BP algorithm 
is also similar to the recent approximate message passing 
(AMP) algorithm in (24) , (26) and this connection may be 
useful for further analysis. 

Numerical Results: Through simulations for femtocell deploy- 
ments, we demonstrate that approximate BP provides good 
performance for sub-band and power allocation to maximize 
utilities of rates, sub-band and power allocation to maximize 
a weighted sum of rates, and beamforming optimization to 
maximize utilities of average rates. Also, although the BP 
algorithm requires multiple exchange of messages before each 
scheduling decision, our simulations indicate good perfor- 
mance with only two rounds of messaging. Thus, the approxi- 
mate BP approach is a promising paradigm for new emerging 
cellular deployments with large interference variations in time 
and space compared to current predominantly macro-only 
deployments. 

A. Previous Work 

Algorithms for IOC in LTE macrocells have been con- 
sidered in a large number of works in both the uplink and 
downlink (27)-(29) . These works are generally based on 
adaptive subband scheduling and fractional frequency reuse 
(FFR) methods (18) , (30) and exploit statistics over large 
numbers of mobiles per macrocell. Interference mitigation 
in femtocells has focussed on similar techniques as well 



as frequency planning, power control |31|, or semi-static 
resource allocation (7), (8), (32). As we will demonstrate 
in the simulations, the BP methods presented here can also 
be used for adaptive subband scheduling as a special case 
of the linear mixing interference model. Also, much of the 
ICIC work has considered slowly varying allocations that don't 
change over few 10s to 100s of milliseconds. Due to the low 
messaging overhead, it is possible that approximate BP can 
also be used for more dynamic interference management in 
femtocell deployments, where there is high variability in load 
and interference from one timeslot to another. 

For scheduling based on more dynamic traffic statistics 
such as queue lengths and head-of-line delays, variants of 
maximum weight scheduling can be used (33), (34). Unfor- 
tunately, computing a maximum weight schedule is generally 
NP-hard, and much work has thus focused on approximate 
algorithms. In addition to the works mentioned above, the 
works (33) and (35) proposed randomized linear complexity 
(but centralized) algorithms, and [36|-[38| present simple 
distributed algorithms for combinatorial interference models. 
Greedy maximal weight matching for such interference models 
has been considered in [39|-]4"T). 



However, many of the above works apply a hard constraint 
interference models, where neighboring links cannot transmit 
simultaneously. Cellular systems in contrast permit multiple 
interfering links to transmit simultaneously and then use rate 
control to adapt to the resulting signal-to-interference and 
noise ratio (SINR). Thus, the degradation in rate with interfer- 
ence is gradual, and are difficult to capture in the combinatorial 
interference model. In contrast, "soft" interference effects can 
be easily modeled in the BP utility framework. 

However, we note that there is an important theoretical 
connection between the methods in [36], (37) and the BP 
method considered here. As we will discuss below, BP arrives 
at the scheduling decision by estimating the marginals of a 
certain joint probability distribution function given in (7). The 
CSMA-type methods in (36) , (37) can be seen as a simulated 
annealing (SA) method for selecting a random scheduling 
vector from precisely the same distribution in the context of a 
constraint combinatorial interference model. SA can be seen 
as an asymptotically exact but slow method (20) for solving 
the optimization problem. In constrast, BP is approximate, but 
potentially faster. 

General overviews of BP can be found in a number of works 
including (42), (43). In the context of wireless scheduling, 
theoretical guarantees have been obtained for on-off channels 
and the combinatorial contention graph model [44], [45]. The 
methods here can be seen as a generalization of these methods 
to soft interference models with larger class of scheduling 
vectors. 

II. Problem Formulation 

We consider a wireless scheduling problem with n links. 
The transmitter of each link j = 1, . . . , n, denoted TX j, must 
select some scheduling vector ■ G X C E" x , which contains 
n x parameters related to link i. Examples of these parameters 
will be given below. The selection of the scheduling vectors 
results in an interference vector Zj £ W lz at the receiver of 
each link i, denoted RX i. The interference is assumed to be 
a linear function of the scheduling vectors of the other links, 



(i) 



for some matrices Ay £ ]g>™z xr ^. We assume that An = 0, 
so that link i does not interfere with itself. We will let x and 
z be the column vectors with entries Xj and z^, 

x = ft ■..<]' eR"-" 
z = [zi.-.z'j'eR"*", 

and write z = Ax where A is the block matrix with entries 
Aij . We call A the interference matrix. Also, we let A,; denote 
the ith column of the A so that z^ = A^x. 

Associated with each link i, is some utility function 
fi(xi,Zi) of the scheduling vector x; and interference vector 
Zi. The scheduling problem is to maximize the overall utility 



maxF(x), -F(x) 



i=i 



/i( x ij z i)- 



(2) 
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We will sometimes call the optimization problem pi, an 
optimization with linear mixing to stress the linear dependence 
of the interference on the transmit vectors. 

III. Linear Mixing Utility Examples 

The linear mixing formulation above is extremely general 
and can incorporate a large class of utility functions and 
interference models. 

A general treatment of utility functions for wireless schedul- 
ing can be found in |46) , (47). In our simulations, we will 
consider utility maximization for both static and time-varying 
problems. For static optimization, the scheduling vectors Xj 
are selected once for a long time period and the utility function 
is typically of the form 



/i(xi,Zi) = Ui(Ri(xi,Zi)), 



(3) 



where Ri(x i: Zi) is the long-term rate as a function of the 
TX vector Xi and interference Zi and Ui(R) is the utility as 
a function of the rate. The problem formulation above can 
incorporate any of the common utility functions including: 
Ui(R) = R which results in a sum rate optimization; Ui(R) = 
\og(R) which is the proportional fair metric and Ui(R) = 
—f)R~P for some (3 > called an /3-fair utility. Penalties can 
also be added if there is a cost associated with the selection 
of the TX vector xj such as power. 

To accommodate time-varying channels and traffic loads, 
many cellular systems enable fast dynamic scheduling in time 
slots in the order of 1 to 2 ms. For these systems, the utility 
maximization can be re-run in each time slot. One common 
approach is that in each time slot t = 0, 1, 2 . . ., the scheduler 
uses a utility of the form 



fi(t,-X.i,Zi) = W i (t)i?j(t,X i ,Z i ), 



(4) 



where Wi(t) is a time-varying weight given by the marginal 
utility 

dUi(Ri(t)) 



OR 



(5) 



and Ri(t) is exponentially weighted average rate updated as 

Ri(t+1) = [l-a)Ri{t) + aRi{t,Xi(t),%(t)), (6) 

where X;(£) is the TX vector and Zj(t) is the interference for 
link i at time t. Any maxima of the optimization |2} with 
the weighted utility Q is called a maximal weight matching. 
A well-known result of stochastic approximation flO) is that 
if a — > 0, and the scheduler performs the maximum weight 
matching with the marginal utilities ((5), then for a large class 
of processes, the resulting average rates will maximize the 
total utility J2i U i(M t ))- 

The above utilities are designed for infinite backlog queues. 
For delay sensitive traffic, one can take the weights Wi(t) to 
be the queue length or head-of-line delay. Maximal weight 
matching performed with these weights generally results in so- 
called throughput optimal performance |9|. These results also 
apply to multihop networks with the so-called backpressure 
weights. 

In addition to incorporating general utilities, an appealing 
feature of the linear mixing framework is that a large class 



of interference models can also be considered, including, for 
example: 

• Flat fading with power control: In this case, Xj is a scalar 
representing the transmit power, and Aij is the gain from 
TX j to RX i, so that Zi is the total interference at RX i. 
The rate, Ri(xi,Zi) can then be described as a function 
of the SINR giXi/zi, where gi is the channel gain along 
link i. Arbitrary SINR to rate mappings may be used. 
Note that a special case of on-off channels where Xj is 
zero or a maximum transmit power can be used. 

• Multiple subbands: The above example is easily extended 
to the case of multiple subbands. As described in the 
Introduction, subband scheduling is one of the key moti- 
vating features of LTE, but the optimization is difficult. 
To handle multiple subbands, we simply let xj and Zi be 
the vectors of transmit and interference powers in each 
subband and Ay be a diagonal matrix with channel gains 
in each subband. 

• Beamforming and linear preceding: The linear mixing 
formulation can also incorporate problems with transmit 
beamforming or linear precoding. For example, suppose 
a link has N transmit antennas and one receive antenna. 
If each transmitter TX j uses a beamforming vector g 
C^, and gij G C N is the channel from TX j to RX i, 
the interference at RX i is given by 

which is linear in the rank one matrices bjb^. Hence, if 



we let x., g 



be the column vector with entries of 



the matrix bjb^, the interference Zi can be represented 
as a linear combination of the vectors Xj. The idea can 
also be generalized to precoding matrices with multiple 
transmit streams. 



IV. Belief Propagation 

A. Standard BP 

We begin by briefly reviewing how we would apply standard 
BP to the optimization |2]). Let u > and define the 
probability distribution 



1 1 " 

K x ) = — exp(uF(x)) = — Y[exp(ufi(xi,Zi)) 



(7) 



where Z is a normalization constant called the partition 
function (it is a function of u). BP can be seen as a method 
to estimate the marginal distributions of the distribution p(x) 
with respect to the variables Xj. From these marginals, one can 
compute the marginal expectations Xj = E(xy). A standard 
result of large deviations |48| is that as u — > oo, under suitable 
conditions, p(x) concentrates around the maxima of F{x) and 

lim x = argmaxi^x). 

u — ^oo x 

So, if we can estimate the marginal expectations of the 
probability distribution (j7]) for large u, we can recover a good 
estimate for the maximization of (j2j. 
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To compute the marginal distributions, BP associates with 
the interference matrix A a bipartite graph G = (V, E) called 
the factor or Tanner graph. The vertices V consists of n 
transmitter nodes associated with the transmitters TX j, and 
n receiver nodes associated with the receivers RX i. There is 
an edge G E if and only if i = j or Ay is non-zero - 
that is TX j has some influence on the interference or signal 
at RX i. We let N IX (i) and N tx (j) be the neighbors sets of 
the nodes RX i and TX j in graph G, respectively. 

With this graph, BP iteratively passes beliefs along the edges 
of the graph that represent estimates of the marginal distribu- 
tions of p(x) with respect to the variables Xj. In the context of 
the wireless scheduling problem, we can interpret the iterations 
as rounds, where computations are first performed at the 
receivers and then at the transmitters. We index the round by t, 
and let Pi^j (t,-) : X n- M denote the belief message from RX 
i to TX j in the receiver half of the round. The reverse belief 
message from TX j to RX i is denoted p.^j(t, •) : X i-» R. 
Pi^j(t,x.j) and pi^jftjXj) denote the values of the beliefs 
at Xj. After some fixed number of rounds, the algorithm is 
stopped and a final scheduling decision, meaning a selection 
of the TX vectors Xj, is made by the transmitters. The steps 
for BP are as follows: 

1) Initialization: Set t = and for all € E, let 
Pi4-j(t,x.j) be some initial distribution on Xj. This dis- 
tribution could be, for example, the uniform distribution 
on the set X. 

2) RX node update: In the RX half of the round, each RX 
i sends a belief message to the transmitters TX j with 
j e N r . x (i) given by 

Pi-,j(t,Xj) = E [exp(ufi(xi, Zj)) | Xj] , (8) 

where z, = A^x as in ([T]l and the expectation is over 
independent x fc ~ pu- k (t,x k ), Vfc e iV rx (i). 

3) TX node update: In the TX half of the round, each TX j 
sends a belief message back to the receivers RX i with 
* G Ntx(j) given by 

Pi^j (t + l,xj) = - Y[ pi^j (t, xj ) , (9) 

where Z is a normalization constant and the product is 
over all I £ N tx (j) with £ ^ i. The iteration number is 
incremented, t = t + 1, and we return to step 2 until a 
sufficient number of rounds have been performed. 

4) Final solution: The final estimate for the marginal dis- 
tribution of Xj is given by 

p j (t + l,x j ) = - 11 p^j(t,Xj). (10) 

ieJV tx (j) 

The scheduling vector can be selected as the maximum 

of this marginal distribution. 
In the case when the graph G has no cycles, it can be shown 
that the pj(t,xf) converges to the true marginal distribution 
of p(x) in with respect to the variable Xj. However, for 
general G, BP is approximate. A complete treatment of BP is 
beyond the scope of this work - the reader is referred to the 
references above. 



Implementation of the above BP algorithm in wireless 
networks is challenging due to two reasons: 

• High computational complexity: The expectation in ([8]) 
requires integration over all the variables x r with r E 
N rx (i) and r ^ j. This computation grows exponentially 
in |iV rx (z)|, which is the number of transmitters interfer- 
ing with RX i. If this set is large, the computation is 
prohibitive. 

• High messaging overhead: Passing the beliefs requires 
unicast messages between each RX and each TX (which 
are neighbors as per graph G) as opposed to single 
broadcast message. Also, in each round t, the messages 
comprise of the beliefs pi^j(t,Xj) and pi^j(t,Xj) for 
all values Xj G X. If X is large, the messaging overhead 
may be significant, and it grows with the number of 
transmitters interfering at a receiver. 

B. Gaussian Approximation 

We first consider the simplification of the RX node up- 
date (|8). We describe the simplification in log domain. Let 

Aj_).j(£, •) and A^j(£, ■) be log likelihood functions, meaning 
any functions such that 

Aj_>j(t,Xj) := - log [pi_^(t, Xj)] + const, (11a) 

&u-j(t,Xj) := - log \pa-j{t, Xj)] + const, (lib) 

where the constants do not depend on Xj (although they 
may depend on t and the indices i and j). Observe that we 
can recover the probabilities from the log likelihoods by the 
relation 

Pi_>j(i,Xj) oc exp[uAi->j(t,Xj)] (12a) 
Pi<_j(i,Xj) oc exp[uAj«_j-(t,Xj)] . (12b) 

We now consider the simplification of the RX node update 
([H) under two cases: when i = j and when i ^ j. We 
begin with the case when i = j. The expectation in (|8]l is 
to be evaluated with Zj given by ([T]i with the variables Xj 
being independent and Xj ~ pi^j{t,Xj). Let x,.<_j(i) and 
Qf<-jW be the mean and l/u times the variance of the 
distribution p^jft, ■). Then, under the simplifying assumption 
that the summation in ([T]) consists of a large number of 
independent terms, we can apply the Central Limit Theorem 
and approximate the distribution of z, with 

z i =M(s ii (t),QUt)H, d3) 

where 

s«(*) = zZ A u X ^W ( 14a ) 

which have the interpretation as a mean and variance of the 
interference at RX i. Applying the Gaussian approximation 
( p"3j ) to the expectation ^ shows that we can write Aj_j,j(-) 
in ( |1 1 a| > as 

Ai-xftxi) w - log [.^(x^^), <%(*))] , (15) 
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where Zio(-) is the partition function 

ZioCx^SiijQy) = yexp[uL i0 (x l ,z i ,s ij ;,Q- J ;)]dz i , (16) 
and 

Ifflfx^z^s^.Q-,) = fi{xi,Zi) - -(zj -s i i)'Q^ s (z i -Sj) 

(17) 

and Q u s is the matrix inverse of Qf^ 

Next consider the case i ^ j. A similar argument as above 
shows that, conditional on x,-, the distribution of z s ; can be 
approximated by the Gaussian 



:^(sy(t),Qy(f)/ U ), 



(18) 



where 



/(*) 



E 



A^Xj^ ; 



(t) 



■ Ayx, 



— Sii (t) ~t" Ay (Xj 



Hi- J 



(*)) 



^ ] Ai,.Q i< _ I ,(t)A i? ,, 
reN rx (i)jtj 

q;,.;/) a„q;; ,,:/:a; ; . 



(19a) 



(19b) 



Then, applying the Gaussian approximation ( 18 i along with 
( |1 lb[ ) to the expectation ^ shows that we can write Ai_y(-) 
in ( |1 1 a| > as 

1 log [Z,(A^(i, ■), % (i), Q£(t))] , (20) 



Aj_y(t,x,,-) 
where 



z i{A(')>*ijiQij) = J ex P [ uL i( x i, z i,Sy, Qy)] ebqctej 



and 



L^x^ZjjSj^Qf -) = L i0 (-Ki,Zi,Sij, Q|.) + A(xj 



(21) 



Then, the RX and TX node update steps of the BP algorithm 
with Gaussian approximation of the interference at the RX are: 

• RX node update: The above equations are used to sim- 
plify the standard BP algorithm as follows. Each RX i 
receives mean and variances Xj<_ j (t) and Qf^j (£) of the 
vectors from the transmitters TX j with j 6 N IX (i). RX 
i also receives the function Ai<_j(i,-) from its serving 
transmitter, TX i. RX i then computes the interference 
means and variances in ( fT4| i and ( [T9| ). Then, for each TX j 
it sends back the log likelihoods Aj_y(i, •) by evaluating 
the log partition functions ( fT5] l and ( |20] i. 

• TX noiie update: It can be verified that converting the 
update |9]) to log domain yields 



Ai+-j(t + l,Xj 



E 



(22) 



Each TX j can first computes the log likelihoods 
Aj^j(i + l,Xj) from the log likelihoods A^_y(t, Xj) 
from the receivers RX I with I £ N tx (j), £ ^ i. Then, 
using the log likelihoods Ai«_,-(i+l, x,-) , TX j computes 
the mean and variance Xj<_j(t+1) and Qf<_,-(i+l) from 
the probability distribution p^j(t+l, x 3 ) in ( |12b| i. Then, 



TX j sends messages to the receivers as described in the 

RX node update above. 
After a sufficient number of rounds, TX j can compute the 
final log likelihood 



E 



(23) 



and then compute the final scheduling vector as the one which 
maximizes the above function. 

We have therefore simplified the standard BP algorithm by 
eliminating the exponential complexity of the RX update |[S), 
and replaced this computation with a Gaussian approximation. 

C. First Order Approximations 

Unfortunately, the Gaussian approximation above does not 
significantly reduce the messaging overhead. Each RX and TX 
must still send separate unicast messages to every TX or RX in 
its neighbor set. Also, although the TX must only send a mean 
and a covariance matrix x^j(i) and Qf^_j(t), the receivers 
must send the entire log likelihood functions A.j_y(i, Xj). 

The messaging overhead can be reduced via selective use of 
first order approximations as follows: Divide the edges G 
E with i 7^ j into two sets - weak and strong - depending 
on whether Ay is small or large. Along a strong edge 
RX i and TX j exchange the full unicast messages described 
above. However, for the weak edges, the messages can be 
replaced with a first order approximation described below. The 
precise classification rule between weak and strong edges is an 
algorithm parameter that can be used to trade off complexity 
and accuracy. 

To describe the first order approximation along the weak 
edges, suppose Ay is small for some edge Let x,(i) 

be the mean value of the distribution corresponding to the log 
likelihood Aj(t, •) in ( |23| ). Now consider the log likelihood 
Ai_^(i,xf) in Applying <[L9) 

A^-(i,x 3 ) « - log [Zi(An-i(t, •),%(*), Qy(*))] 



H 1 



log [ZiiA^iit,- 



u 



■ Aijte -%{t)),QZ&))\ 



(*>) 



,(*) X 3 



const 



(24) 



where (a) follows from ( fT9| ) and the approximation that 
Qy-(t) « Q|j(i) and Xj<_j(t) « Xj(i) when Ay is small; 
and (b) follows from taking a Taylor's approximation with 



uy (t) = A^-D^i) - A^A 2 (i)Ayx,(i), 
where for r = 1,2, D ir (t), are the derivatives 
1 <9 r 



(25) 



.(<) := -7^ log [Zi(A, 



i(-),^i(*),Q'i(*))]. (26) 



The constant in p4| i is independent of Xj. Using standard 
properties of the cumulant function J43) , it can be shown that 
the derivative is given by 



Ar(*) = E 



5 r 



(27) 
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where the expectation is with respect to the conditional distri- 
bution 

ftt lSi |s,Q«(XtjZ< I s„ (*),(%(*)) 

:= — exp[uL i (x i ,z i ,s ii (t),Q| i (f)] , (28) 
A 

and £,(•) is defined in ( |2T| ) and (17) . The dependence on 
in ([27]) is implicit in ( |2"T| ). Note that the derivatives 
Di r (t) in ( |27| ) can be interpreted as a sensitivity of the 
expected utility /j(xj,Zj) to changes in the interference z$. 

The computation of Xj«_j(i) can also be simplified as 
follows: Recall that x.i^j(t) and x,-(i) are the expectation 
with respect to the likelihood functions A^^t, •) in ( |22"] > and 
Aj(<, •) in ( f2"3~| l, respectively. Therefore, we can write them as 

xV;(t) = E [ Xj - 1 A^ (*,•)] (29a) 
Stj(t) - E[x 3 |A,(t,-)], (29b) 

where we have used the notation E(g(x)|A(x)) to denote the 
expectation of g(x) with respect to the probability distribution 

Pa(x) oc exp [tiA(x)] . 

It is easy to check that for small perturbations of e(x) of A(x), 
we have the first order approximation 

E(g(x)|A + e)wE(5(x)|A) 

+ u [E( 5 (x)e(x)|A) - E( ff (x)|A)E(e(x)|A)] . (30) 

Therefore, 

Xi<-i(i + 1) = E [xjlAjit + 1, •) - A^(i, •)] 

(b) 

« E ( Xj - 1 A, (t + 1 , • ) ) - «E ( Xj - Ai_>j (t , x, ) | Aj (t + 1 , • ) ) 
- uE(x,|A,(t, ■))E(A l ^-(t,x,)|A,(t + 1, •)) 

« £,-(« + l)-Qf(t + l)uy(t) (31) 

where (a) follows from ( |29a| i along with ( p2) and ( |23| ); (b) 
follows from pO) ; and (c) follows from ( |29b| i and ( |24] >. 

We can use the above relations to define the following 
algorithm: 

1) Initialization: Set t — 0. Each TX j broadcasts an initial 
scheduling vector Xj-(t) and variance QJ(t). These can 
be based on the mean and variance of Xj over the set X. 
The receivers RX i sets Xi«-j(t) = Xj(i) and Qf<_j(i) = 
Q|(i) for all j, and Aj<_j(t,Xj) = for all x 4 . 

2) /?X nocfe update: In the RX half of the round, each RX 
i first computes the interference means and variances 
sa(t) and Qf i (t) in ([14]) and sy(i) and Qy(i) in 
( fl9) , Then, RX i computes the log likelihood function 
Aj_j,j(t, x,-) in (JT3J and sends it as a unicast message 
to its serving transmitter TX i. Also, for each strong 
edge RX i computes the log likelihood functions 
Ai->,j(t, Xj) in ( |20] l and sends it as a unicast message to 
the interfering TX j. For the weak edges, RX i simply 
computes the sensitivity Dn(t) in ( |27] i and broadcasts 
it to all other transmitters TX j. The receiver also 
computes (t) in ( |25] > and stores it for the next round. 
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Fig. 1. One round of the BP messages along the weak edges interpreted as 
a "soft" RTS / CTS mechanism. 



3) TX node update: In the TX half of the round, each 
transmitter TX j would have received the log likelihoods 
Aj_j,j(t,Xj) from the receivers RX i, along the edges 

that were strong. For any weak edge TX j 

can approximately compute the log likelihood from the 
sensitivity Dn(t) using ( p4] i. TX j can then compute 
the log likelihoods Ai^-j(t + l,Xj) from ( f22] ) for all 
i G N t x (j) and tne l°g likelihood Aj(t + l,Xj) in ( |23] i. 
TX j sends the receiver RX j that it is serving the entire 
log likelihood Aj<_j(t + l,Xj). For the receivers RX i 
such that is a strong edge, TX j computes the 

mean Xj«_j-(i + 1) and variance Qf^j(t + 1) from the 
log likelihood A,- < _ J -(t+ l,Xj) and sends it as a unicast 
message to RX i. Each TX j also computes the mean and 
variance x.j(t + l) and QJ(i+l) from the log likelihood 
Aj(t+ l,Xj) and broadcasts the quantities to the other 
receivers. Any receiver RX i that is along a weak edge 

can then approximately compute Xi^-j(t) from 

The round number is incremented, £ = t + 1, and we 
return to Step 2 until a fixed number of rounds have 
been performed. 

4) F/na/ solution: After the final round, each transmitter 
TX j takes the scheduling vectors to be the vector Xj 
that maximizes the log likelihood Aj(t, X,-). 

The message flow along the weak edges has an appeal- 
ing interpretation. Consider Fig. [T] where a transmitter TX1 
attempts to send data to a receiver RX1. The receiver RX1 
experiences interference from a second transmitter TX2, while 
the transmitter TX1 causes interference onto a victim receiver 
RX3. Fig. [T] shows the messages along the weak edges in 
one round of the BP algorithm to coordinate the interference. 
The transmitters TX1 and TX2 will broadcast the mean and 
variance Xj(i) and QJ(<) of their intended transmit vectors. 
These transmissions can be interpreted as "soft" request to 
sends (RTS). They are soft since the intended transmit vectors 
are signaled by a distribution. Based on the transmit vector 
distribution from the interfering TX2, the receiver RX1 replies 
to the serving TX1 with an estimate of the interference s 1:L (t) 
and Qfi(i) which can be interpreted as a "soft" channel 
quality indication (CQI). As victim receivers, RX1 and RX3 
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Method 


RX i 


TX j 


Exact 


0(|Af|l JV »»l) 


0(\X\\N tx U)\) 


Gaussian Approx. 




0(\X\\N tK (j)\) 


First Order 


0{\Nr*{i)\) 


0(|iVtx(i)|) 




TABLE I 




Computational Complexity 


Method 


RX i 


TX j 


Exact 


0{\X\\N^(i)\) 


0(|^||JV tx (j)|) 


Gaussian Approx. 


0(\X\\N IX (i)\) 


o(|JVteO')l) 


First Order 


0(1) 


O(l) 



TABLE II 
Communication Complexity 



also compute the sensitivities to the interference level by the 
derivatives Di r (t). These values can be interpreted as soft clear 
to send (CTS) indications, since instead of a binary go/no go 
type CTS, they signal a soft cost on changes in the interference 
from other transmitters. 

D. Communication and Computation Complexity 

We summarize the complexity of one round of the three 
variations of BP in Tables [j] and [TT] Specifically, we focus on 
RX i and TX j and consider how the complexity grows with 
X, JV«(i) and AT tx (j). 

V. Numerical Simulation 

The BP algorithm was simulated on a a simplified version 
of an industry standard model for LTE femtocell evaluation 
in j3). The simulation parameters are shown in Table III The 



network consists of a 3 x 3 grid of 10m x 10m apartments 
with active links in 5 of the 9 apartments. Each link consists of 
one femto BS transmitting to one femto mobile (called a UE, 
or user equipment, in 3GPP terminology). Due to restricted 
association, UEs connect to the femto BS in their apartment 
even if it is not the BS with the minimum path loss. As 
mentioned in the introduction, this scenario exposes many 
links to strong interference, and thus presents a good test 
scenario for advanced interference coordination algorithms. 

In the first simulation, we considered a time-varying sim- 
ulation with a simple on-off model where, in each time slot, 



Parameter 
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Network topology 


3x3 apartment model, with active links in 
5 of the 9 apartments. 


Carrier freq 


2 GHz 


Bandwidth 


5 MHz 


Wall loss 


or 10 dB 


Lognormal shadowing 


10 dB std. dev. 


Path loss 


38.46 + 20 log 10 (_R) + 0.7R dB, R distance 
in meters. 


Femto BS TX power 


dBm 


Femto UE noise figure 


4 dB 
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Fig. 2. Downlink femtocell simulation with an on-off channel model and 
time-varying flat fading. The top panel shows the rate CDFs across all links 
and drops using various optimization algorithms for weight matching. The 
bottom panel is the CDF of system utility (represented as the harmonic mean 
rate) across different drops. 



each link either transmits at the max power or is completely 



off. As described in Section III for the time-varying problem 



the utility maximization was rerun in each time slot with 
the weighted sum rate utility Q with weights (15). We used 
the proportional fair utility Ui(R) = \og(R). We generated 
independent flat fades on the links in each slot, and took a filter 
time constant in |6]i was a = 0.1. For each random realization, 
or "drop", of the femto network, we ran the simulation over 
100 time slots and measured the average rate. The wall loss 
in the femto model was 10 dB. 

The top panel of Fig. [2] plots the cumulative distribution 
function (CDF) of the time-averaged rates for 5 links and 100 
drops comparing various optimization methods for computing 
the maximum weighted matching optimization (p). The curve 
"reuse 1" is the case when all links transmit at max power. 
Two cases of BP are simulated: (i) 4 rounds of BP in each 



TABLE III 
Simulation parameters. 



time slot using the Gaussian approximation in Section IV-B but 
no first order approximations; and (ii) using only two rounds 
of BP with both the Gaussian approximations and first order 
approximations on all interfering links less than dB below 
the serving link. We see that even the linear approximated BP 
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Fig. 3. Subband static optimization. Performance comparison for downlink 
femtocell rates with various optimizations for 4 independently faded subbands. 



with 2 rounds does significantly better than reuse 1, and is 
not that far from the 4 round BP. The gap between BP and 
reuse 1 particularly large at low rates. For example, for the 
20% worst links, BP offers almost a factor of 5 improvement 
in rate over reuse 1 . 

Also plotted is the rate CDF with optimal matching in 
each time slot based on an exhaustive centralized search. BP 
performs reasonably close to this curve, although there is still 
an obvious gap, again at low rates. 

The bottom panel of Fig. [2] shows the CDF of the total 
system utility over the 100 drops. Since the optimization used 
a log utility, the optimization is equivalent to maximizing the 
harmonic mean rate over the 5 links. We see in this plot 
that the approximate BP algorithms are close to optimal and 
significantly better than reuse 1. 

As a second simulation, we considered the same problem 
but with a static single optimization and K = 4 subbands. 
Random fading was generated in each subband, so the simu- 
lation would be applicable in the case with frequency selective 
fading with coherence bandwidth roughly equal to the subband 
bandwidth. The optimization was over scheduling vectors 



x, e 



with each component being on or off, so there are 
2 K — 1 = 15 possible non-zero scheduling vectors in each 
link. Optimal subband allocation is a well-known challenging 
but important problem for OFDMA systems like LTE. 

Under these assumptions, Fig. [3] shows the CDF of the rates 
under various optimization methods. Similar to the dynamic 
single subband case, we see that BP provides significant gains 
over simple reuse 1, even when we use linear approximations 
and only four rounds. Also, for most links, BP achieves a rate 
reasonably close to the optimal subband allocation found by 
exhaustive search over all subband allocations over all rates. 
However, for the lowest rate links, there is a significant gap 
between BP and optimal. Thus, even though BP outperforms 
reuse 1 significantly in this regime, there is significant room 
for improvement. 

The BP method can also be applied to beamforming (BF) 
problems. As a simple simulation, we consider again the 



Fig. 4. Beamforming optimization. Downlink femtocell rates for various 
optimizations with 2 TX antenna beamforming. Channel models assume 
no scattering and beamforming optimization is performed over linear phase 
beamforming vectors. 



femtocell deployment with transmit beamforming with two 
antennas with half wavelength spacing spacing and one receive 
antenna. We neglect scattering so the channel appears as a 
linear phase across the two antennas on all links. For the 
transmit vectors Xj, we optimize beamforming angles over 
10 angles uniformly spaced between and tt. A performance 
comparison of the rate CDFs under various optimization 
algorithms is shown in Fig. |4] In this case, we took a wall loss 
of dB. The curve labeled "opt serving link only" is the case 
when the BF vector is chosen to maximize the signal strength 
from the serving link only without regard to interference. This 
simple method provides the baseline. We see that BP provides 
some gains over serving link optimization. For example, the 
median rate with BP is approximately 50% higher than using 
optimal BF on the serving link only. Also, BP appears to 
be reasonably close to the optimal BF selection based on 
exhaustive search. 

VI. Conclusions 

We formulated a general wireless scheduling and interfer- 
ence coordination problem as an optimization problem with 
linear mixing utilities. This was cast in a BP framework 
where the goal is to compute the marginal distributions of a 
joint probability function. Using Gaussian and linear approx- 
imations, we obtained a distributed interference coordination 
algorithms with low overhead. The algorithm has a natural in- 
terpretation as a soft RTS/CTS scheme. Numerical simulations 
demonstrated that the resulting algorithm is close to an optimal 
scheme for dynamic and static sub-band optimization as well 
as for beamforming coordination across cells. Moreover, the 
results show that the algorithm computes a good operating 
point in just two to four iterations making it very attractive to 
be used in practical wireless cellular systems. In the future, 
we plan to explore connections between our work and the 



AMP framework in \22\ to obtain performance bounds for 
large random networks. 
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